{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GIAO 的 RHF 核磁 (屏蔽) 共振常数 (NMR) 数值导数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 创建时间：2020-08-31"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这篇文档中，我们会讨论使用 PySCF 以及其作为 libcint 的接口，计算 GIAO 的 RHF 数值核磁 (屏蔽) 共振常数 (Nuclear Magnetic Resonance constant, NMR) 的程序。该文档大量参考 PySCF 的代码 [nmr/rhf.py](https://github.com/pyscf/pyscf/blob/master/pyscf/prop/nmr/rhf.py)。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "与上一篇文档一样，我们的讨论中所使用到的分子体系 `mol` 会是非对称的氨分子，并且取用最小基组。其 RHF 计算放在实例 `mf`，而 NMR 计算实例会放在 `mf_nmr`。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscf import gto, scf, dft\n",
    "from pyscf.prop import nmr\n",
    "from pyscf.data import nist\n",
    "from scipy import constants\n",
    "import numpy as np\n",
    "\n",
    "np.set_printoptions(precision=5, linewidth=150, suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol = gto.Mole()\n",
    "mol.atom = \"\"\"\n",
    "N  0.  0.  0.\n",
    "H  0.  1.  0.2\n",
    "H  0.1 0.3 1.5\n",
    "H  0.9 0.4 -.2\n",
    "\"\"\"\n",
    "mol.basis = \"STO-3G\"\n",
    "mol.verbose = 0\n",
    "mol.build()\n",
    "coord_orig = np.zeros(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "nocc, nao, nmo, natm = mol.nelec[0], mol.nao, mol.nao, mol.natm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其自洽场能量为"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-55.25354051468657"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf = scf.RHF(mol).run()\n",
    "mf.e_tot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "核磁屏蔽张量 (Shielding Constant) $\\sigma_{ts}^A$ 可以表示如下 (维度 $(A, t, s)$)：\n",
    "\n",
    "$$\n",
    "\\sigma_{ts}^A = \\frac{\\partial^2 E_\\mathrm{tot}}{\\partial \\mathscr{B}_t \\partial \\mu_{A_s}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[305.25641,  22.67382,  12.29714],\n",
       "        [-22.38193, 180.11491,  -2.84064],\n",
       "        [ 13.65836, -28.07439, 383.02872]],\n",
       "\n",
       "       [[ 25.68677,  -0.27932,  -0.31477],\n",
       "        [ -0.84315,  35.70697,  -4.26016],\n",
       "        [ -1.06334,   2.07184,  28.17535]],\n",
       "\n",
       "       [[ 24.06192,  -0.12066,   0.81476],\n",
       "        [  0.99821,  24.76595,  -0.86457],\n",
       "        [ -0.74954,   0.2574 ,  31.65924]],\n",
       "\n",
       "       [[ 39.48533,   4.28363,  -7.95068],\n",
       "        [  0.89957,  26.69818,  -0.9659 ],\n",
       "        [ -4.40617,  -0.48147,  28.90276]]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mf_nmr = nmr.RHF(mf)\n",
    "mf_nmr.kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上述结果是以 ppm 为单位的表达。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 微扰原子矩阵的表达"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 外磁场微扰的一阶 Core Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "关于该微扰矩阵，我们已经在前两篇文档中有比较多的说明了。\n",
    "\n",
    "$$\n",
    "\\begin{split}\\begin{align}\n",
    "h_{\\mu \\nu}^{\\mathscr{B}_t}\n",
    "&=\n",
    "\\frac{1}{2} \\langle \\mu | \\hat l_t | \\nu \\rangle_{\\mathrm{Gauge} \\rightarrow \\boldsymbol{R}_\\nu}\n",
    "+ \\langle U_\\mathrm{g}^t \\mu | \\hat t | \\nu \\rangle\n",
    "+ \\langle U_\\mathrm{g}^t \\mu | \\hat v_\\mathrm{nuc} | \\nu \\rangle\n",
    "\\\\ & \\quad\n",
    "+ \\sum_{\\kappa \\lambda} ( U_\\mathrm{g}^t \\mu \\nu | \\kappa \\lambda ) D_{\\kappa \\lambda}^{(0)}\n",
    "- \\frac{1}{2} \\sum_{\\kappa \\lambda} ( U_\\mathrm{g}^t \\mu \\lambda | \\kappa \\nu ) D_{\\kappa \\lambda}^{(0)}\n",
    "- \\frac{1}{2} \\sum_{\\kappa \\lambda} ( U_\\mathrm{g}^t \\kappa \\nu | \\mu \\lambda ) D_{\\kappa \\lambda}^{(0)}\n",
    "\\end{align}\\end{split}\n",
    "$$\n",
    "\n",
    "由于在 PySCF 中，函数 `nmr.rhf.make_h10` 就是专门用于生成该矩阵的，因此我们就用下述代码生成 `hcore_1_B` $h_{\\mu \\nu}^{\\mathscr{B}_t}$ 表示该矩阵。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "dm_guess = mf.make_rdm1()\n",
    "hcore_1_B = 1j * nmr.rhf.make_h10(mol, dm_guess)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 外磁场微扰的一阶重叠矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "该微扰矩阵也已经在前两篇文档中有比较多的说明了。\n",
    "\n",
    "$$\n",
    "S_{\\mu \\nu}^{\\mathscr{B}_t} = \\langle U_g^t \\mu | \\nu \\rangle\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们用下述代码生成 `ovlp_1_B` $S_{\\mu \\nu}^{\\mathscr{B}_t}$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ovlp_1_B = - 1j * mol.intor(\"int1e_igovlp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 核磁偶极的一阶 Core Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "核磁偶极 $\\mu_{A_t}$ 所产生的一阶算符贡献可以表达为\n",
    "\n",
    "$$\n",
    "\\hat h {}^{(1)} (\\boldsymbol{\\mu}_A)\n",
    "= - i \\alpha^2 \\boldsymbol{\\mu}_A \\cdot \\left( \\boldsymbol{\\nabla} \\frac{1}{\\boldsymbol{r - \\boldsymbol{R}_A}} \\times \\boldsymbol{\\nabla} \\right)\n",
    "= - i \\alpha^2 \\boldsymbol{\\mu}_A \\cdot \\left( \\frac{\\boldsymbol{r} - \\boldsymbol{R}_A}{|\\boldsymbol{r} - \\boldsymbol{R}_A|^3} \\times \\boldsymbol{\\nabla} \\right)\n",
    "$$\n",
    "\n",
    "其中，$\\boldsymbol{R}_A$ 表示原子 $A$ 的核坐标，$\\alpha$ 表示精细结构常数，$1/\\alpha \\simeq 137$。该常数可以从 PySCF 中获得，也可以从 SciPy 中获得。注意等式左边的 $\\boldsymbol{\\mu_A}$ 看作是外加的核磁偶极大小，而等号右边的 $\\mu$ 表示原子轨道，两者意义不同；等号左边角标 $\\boldsymbol{A}$ 表示原子核坐标向量，而之前两篇文档中的 $A$ 在很多文章或教材中表示 $\\frac{1}{2} \\boldsymbol{B} \\times \\boldsymbol{r}$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "137.0359990836958"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 / constants.alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "但为了方便，在最后核算结果之前，我们会暂且将 $\\alpha$ 当作 1 来处理。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在 PySCF 中，实现上述过程的积分字符是 `int1e_ia01p`；但其使用需要告知 `gto.mole.intor` 函数以其 $1 / \\boldsymbol{r}$ 的规范原点位置具体处在哪个原子核中心。\n",
    "\n",
    "$$\n",
    "h_{\\mu \\nu}^{\\mu_{A_t}} = - i \\langle \\mu | \\boldsymbol{\\nabla} \\frac{1}{\\boldsymbol{r}} \\times \\boldsymbol{\\nabla} | \\nu \\rangle_{\\text{Gauge of } \\boldsymbol{r} \\rightarrow \\boldsymbol{R}_A}\n",
    "$$\n",
    "\n",
    "我们用下述代码生成 `hcore_1_m` $h_{\\mu \\nu}^{\\mu_{A_t}}$：(维度为 $(A, t, \\mu, \\nu)$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "hcore_1_m = np.zeros((natm, 3, nao, nao), dtype=np.complex128)\n",
    "for atom_idx in range(natm):\n",
    "    with mol.with_rinv_orig(mol.atom_coord(atom_idx)):\n",
    "        hcore_1_m[atom_idx] = - 1j * mol.intor(\"int1e_ia01p\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 磁场与核磁偶极的二阶 Core Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "磁场与核磁偶极之间的算符乘积会产生二阶算符贡献项：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat h {}^{(2)} (\\boldsymbol{\\mathscr{B}}, \\boldsymbol{\\mu}_A) | \\nu \\rangle\n",
    "&= \\frac{\\alpha^2}{2} \\boldsymbol{\\mathscr{B}}^\\mathrm{T} \\left( (\\boldsymbol{r} - \\boldsymbol{R}_\\nu) \\cdot \\boldsymbol{\\nabla} \\frac{1}{\\boldsymbol{r} - \\boldsymbol{R}_A} - (\\boldsymbol{r} - \\boldsymbol{R}_\\nu) \\boldsymbol{\\nabla} \\frac{1}{(\\boldsymbol{r} - \\boldsymbol{R}_A)^\\mathrm{T}} \\right) \\boldsymbol{\\mu}_A | \\nu \\rangle \\\\\n",
    "&\\quad + \\alpha^2 \\boldsymbol{\\mathscr{B}}^\\mathrm{T} \\boldsymbol{U}_\\mathrm{g} \\left( \\boldsymbol{\\nabla} \\frac{1}{\\boldsymbol{r} - \\boldsymbol{R}_A} \\times \\boldsymbol{\\hat p} \\right)^\\mathrm{T} \\boldsymbol{\\mu}_A | \\nu \\rangle\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "去除精细结构常数 $\\alpha$ 的贡献后，其矩阵的表达形式则是\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "h_{\\mu \\nu}^{A_s t} \\mathscr{B}_t \\mu_{A_s}\n",
    "&= \\langle \\mu | - \\frac{1}{2} \\frac{(t - t_\\nu) (s - s_A)}{|\\boldsymbol{r} - \\boldsymbol{R}_A|^3} | \\nu \\rangle\n",
    "- \\delta_{ts} \\langle \\mu | - \\frac{1}{2} \\sum_{w} \\frac{(w - w_\\nu) (w - w_A)}{|\\boldsymbol{r} - \\boldsymbol{R}_A|^3} | \\nu \\rangle \\\\\n",
    "&\\quad + \\langle U_\\mathrm{g}^t \\mu | \\left( \\boldsymbol{\\nabla} \\frac{1}{\\boldsymbol{r} - \\boldsymbol{R}_A} \\times \\boldsymbol{\\hat p} \\right)_s | \\nu \\rangle\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们用下述代码生成 `hcore_2` $h_{\\mu \\nu}^{A_s t}$：(维度为 $(A, t, s, \\mu, \\nu)$，注意维度 $t$ 对应外磁场，而 $A, s$ 对应核磁偶极)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "hcore_2 = np.zeros((natm, 3, 3, nao, nao))\n",
    "for atom_idx in range(natm):\n",
    "    with mol.with_rinv_origin(mol.atom_coord(atom_idx)):\n",
    "        hcore_2[atom_idx] += mol.intor(\"int1e_giao_a11part\").reshape((3, 3, nao, nao))\n",
    "        hcore_2[atom_idx] -= np.einsum(\"ts, uv -> tsuv\", np.eye(3), mol.intor(\"int1e_giao_a11part\").reshape((3, 3, nao, nao)).trace(axis1=0, axis2=1))\n",
    "        hcore_2[atom_idx] += mol.intor(\"int1e_a01gp\").reshape((3, 3, nao, nao))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数值导数求 NMR 核磁屏蔽张量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "随后我们就可以通过数值梯度求核磁屏蔽张量 $\\sigma_{ts}^A$ (维度 $(A, t, s)$)：\n",
    "\n",
    "$$\n",
    "\\sigma_{ts}^A = \\frac{\\partial^2 E_\\mathrm{tot}}{\\partial \\mathscr{B}_t \\partial \\mu_{A_s}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在此之前，我们仍然需要构造一个通过更改 `get_hcore` Core Hamiltonian 与 `get_ovlp` 重叠矩阵的 PySCF 自洽场实例，以施加外场获得能量的函数 `eng_nmr_field`。其输入的参数 `dev_xyz_B` 是三维外加磁场大小 (对应维度 $t$)，`dev_xyz_m` 是三维外加核磁偶极大小 (对应维度 $s$)，`atom_idx` 是原子序号 (对应维度 $A$)。\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "h_{\\mu \\nu} &= h_{\\mu \\nu}^{(0)} + \\mathscr{B}_t h_{\\mu \\nu}^{\\mathscr{B}_t} + \\mu_{A_s} h_{\\mu \\nu}^{\\mu_{A_s}} + \\mathscr{B}_t \\mu_{A_s} h_{\\mu \\nu}^{\\mathscr{B}_t \\mu_{A_s}} \\\\\n",
    "S_{\\mu \\nu} &= S_{\\mu \\nu}^{(0)} + \\mathscr{B}_t S_{\\mu \\nu}^{\\mathscr{B}_t}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "dm_guess = mf.make_rdm1()\n",
    "\n",
    "def eng_nmr_field(dev_xyz_B, dev_xyz_m, atom_idx):\n",
    "    mf = scf.RHF(mol)\n",
    "    def get_hcore(mol_=mol):\n",
    "        hcore_total  = np.asarray(scf.rhf.get_hcore(mol_), dtype=np.complex128)\n",
    "        hcore_total += np.einsum(\"tuv, t -> uv\", hcore_1_B, dev_xyz_B)\n",
    "        hcore_total += np.einsum(\"tuv, t -> uv\", hcore_1_m[atom_idx], dev_xyz_m)\n",
    "        hcore_total += np.einsum(\"tsuv, t, s -> uv\", hcore_2[atom_idx], dev_xyz_B, dev_xyz_m)\n",
    "        return hcore_total\n",
    "    def get_ovlp(mol_):\n",
    "        ovlp_total  = np.asarray(scf.rhf.get_ovlp(mol_), dtype=np.complex128)\n",
    "        ovlp_total += np.einsum(\"tuv, t -> uv\", ovlp_1_B, dev_xyz_B)\n",
    "        return ovlp_total\n",
    "    mf.get_hcore = get_hcore\n",
    "    mf.get_ovlp  = get_ovlp\n",
    "    return mf.kernel(dm=dm_guess)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "随后使用与前文类似的数值梯度方式就能得到核磁屏蔽常数了；所使用的数值差分大小对于外磁场与核磁偶极均为 $10^{-4}$ a.u.。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "interval = 1e-4\n",
    "num_nmr = np.zeros((natm, 3, 3))\n",
    "for atom_idx in range(natm):\n",
    "    for t in range(3):\n",
    "        for s in range(3):\n",
    "            dev_xyzs_B, dev_xyzs_m = np.zeros((2, 3)), np.zeros((2, 3))\n",
    "            dev_xyzs_B[0, t] = dev_xyzs_m[0, s] = -interval\n",
    "            dev_xyzs_B[1, t] = dev_xyzs_m[1, s] =  interval\n",
    "            num_nmr[atom_idx, t, s] = (\n",
    "                + eng_nmr_field(dev_xyzs_B[0], dev_xyzs_m[0], atom_idx)\n",
    "                - eng_nmr_field(dev_xyzs_B[1], dev_xyzs_m[0], atom_idx)\n",
    "                - eng_nmr_field(dev_xyzs_B[0], dev_xyzs_m[1], atom_idx)\n",
    "                + eng_nmr_field(dev_xyzs_B[1], dev_xyzs_m[1], atom_idx)\n",
    "            ) / (4 * interval**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 5.73238,  0.42579,  0.23092],\n",
       "        [-0.4203 ,  3.38236, -0.05334],\n",
       "        [ 0.25651, -0.52722,  7.19284]],\n",
       "\n",
       "       [[ 0.48237, -0.00525, -0.00591],\n",
       "        [-0.01584,  0.67054, -0.08   ],\n",
       "        [-0.01997,  0.03891,  0.5291 ]],\n",
       "\n",
       "       [[ 0.45185, -0.00226,  0.01529],\n",
       "        [ 0.01875,  0.46508, -0.01624],\n",
       "        [-0.01408,  0.00484,  0.59453]],\n",
       "\n",
       "       [[ 0.74149,  0.08044, -0.14931],\n",
       "        [ 0.01689,  0.50136, -0.01814],\n",
       "        [-0.08274, -0.00904,  0.54276]]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_nmr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "留意到我们之前一直都使用去除结构精细常数的结果，因此我们需要乘以 $\\alpha^2$。同时由于单位是 ppm，因此最终我们需要乘以 $10^6 \\alpha^2$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[305.25688,  22.67379,  12.2969 ],\n",
       "        [-22.38171, 180.11542,  -2.84064],\n",
       "        [ 13.65948, -28.0751 , 383.02872]],\n",
       "\n",
       "       [[ 25.68659,  -0.27935,  -0.3146 ],\n",
       "        [ -0.84328,  35.70701,  -4.26026],\n",
       "        [ -1.06338,   2.07184,  28.17533]],\n",
       "\n",
       "       [[ 24.06171,  -0.12057,   0.81443],\n",
       "        [  0.99832,  24.76594,  -0.86458],\n",
       "        [ -0.74969,   0.25765,  31.65949]],\n",
       "\n",
       "       [[ 39.48546,   4.28353,  -7.95088],\n",
       "        [  0.89964,  26.69806,  -0.9661 ],\n",
       "        [ -4.40616,  -0.48158,  28.9029 ]]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_nmr * constants.alpha**2 * 10**6"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
